% I've (fede) modified this function and eliminated redundancies and
% adapted to Aiyagari files

function  [mu_t0] = bc_update_mut0(r, W_cc, kappa, y_t, a_t, sigma_c,psi, cov_fun, ex_signals, eps_eta,b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c)

%  inputs
%  parameters: r, ybar, rho_y, W_cc, kappa
%  y_t = current state y_t
%  a_t = current state a_t
%  y_hist = past states y_1 ~ y_(t-1)
%  a_hist = past states a_1 ~ a_(t-1)
%  eta_hist = past signals eta_1 ~ eta_(t-1)
%  sigma_etasq_hist = past signal precision^(-1), after forgetting
%  prior functions and parameters: mu0, sigma_c, psi, cov_fun, ex_signals, res_action

% NOTE: mu0 must be the time-zero prior policy function as well as the
% "correct" policy function

%  1) Beliefs based on eta_{t-1}

%  Ignore any point with no information
%y_hist = y_hist(sigma_etasq_hist < 1e10);
%a_hist = a_hist(sigma_etasq_hist < 1e10);
%eta_hist = eta_hist(sigma_etasq_hist < 1e10);
%sigma_etasq_hist = sigma_etasq_hist(sigma_etasq_hist < 1e10);



mu0_ytat = fast_interp( y_t + (1+r)*a_t, ctilde, stateGrid_step, stateGrid,stateGrid_size);

%  prior on chat(y) and Sigma(y,y')
if isempty(hist_mat)
    
    Sigma_t0 = sigma_c^2; %Var(cstar(y) | eta_hist )
    mu_t0 = mu0_ytat;
    
else
    
    y_hist   = hist_mat(:,1);
    a_hist   = hist_mat(:,2);
    eta_hist = hist_mat(:,3);
    sigma_etasq_hist = hist_mat(:,4);
    
    if cov_fun == 1
        
        temp1_yat =  (1+r)*a_t + y_t - (1+r)*a_hist' - y_hist';
        Sigma12  = sigma_c^2*exp(-psi*(temp1_yat.^2));
        
        
        v = Linv*Sigma12';
        sum_vsq = sum(v.^2);
        Sigma_t0 = sigma_c^2 - sum_vsq;
        
        mu_t0 = mu0_ytat + Sigma12*mu_t0_hist;
        
    else
        temp1_yat = bsxfun(@minus, (1+r)*a_t + y_t, (1+r)*a_hist' + y_hist');
        Sigma12  = sigma_c^2*exp(-psi*(abs(temp1_yat)));
        
        
        v = Linv*Sigma12';
        sum_vsq = sum(v.^2);
        Sigma_t0 = sigma_c^2 - sum_vsq;   %Var(cstar(y) | eta_hist )
        
        mu_t0 = mu0_ytat + Sigma12*mu_t0_hist;
    end
    
end



end